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Abstract. The Einstein equations have proven surprisingly difficult to solve 
numerically. A standard diagnostic of the problems which plague the field is 
the failure of computational schemes to satisfy the constraints, which are known 
to be mathematically conserved by the evolution equations. We describe a new 
approach to rewriting the constraints as first-order evolution equations, thereby 
guaranteeing that they are satisfied to a chosen accuracy by any discretization 
scheme. This introduces a set of four subsidiary constraints which are far 
simpler than the standard constraint equations, and which should be more easily 
conserved in computational applications. We explore the manner in which the 
momentum constraints are already incorporated in several existing formulations of 
the Einstein equations, and demonstrate the ease with which our new constraint- 
conserving approach can be incorporated into these schemes. 



PACS numbers: 04.20.-q, 04.25.Dm 

1. Introduction 

Despite significant recent advances in both computational power and algorithmic com- 
plexity, there remain significant unresolved problems with numerical implementations 
of the Einstein equations. Perhaps the most exciting recent developments are the 
many new three-plus-one dimensional formulations of these equations, which, at least 
in part, provide greater stability than the original ADM formulation (see P for a 
review) . 

In particular, formulations of the Einstein equations in strongly hyperbolic, flux- 
conservative form have opened the way for the application of algorithms and techniques 
originally developed for computational fluid dynamics (see, for example, 4 ). Despite 
the great promise and significant advances, modern numerical relativity codes are still 
unable to fully simulate the complete coalescence of binary black hole systems. A 
full understanding and complete simulation of such systems are of vital importance to 
the analysis of gravitational wave signals collected by LIGO and similar gravitational 
wave detectors. 
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Until recently, almost all three-plus-one dimensional formulations of the Einstein 
equations have shared one common feature — the constraint equations are monitored 
through evolution, but never again solved once the initial data is constructed. An 
exception to this is the Texas group, who resolve the elliptic constraint equations 
after every few timesteps The four constraints are more typically treated as 
additional conditions which are monitored while the spacetime is evolved. In fact 
the constraints, and the Hamiltonian constraint in particular, have been found to 
be excellent prognosticators for the accuracy and stability of the numerical solution. 
The rule of thumb appears to be that when the Hamiltonian constraint "explodes" , 
the code will crash shortly thereafter. To the best of our knowledge, all three-plus- 
one dimensional formulations of the Einstein equations suffer from these fundamental 
problems. 

In this paper we propose an approach to transforming all four constraint equations 
into evolution equations for a new set of dynamical variables — essentially the 
conformal factor and it's spatial derivatives, quantities which reduce to the Newtonian 
potential and force in the weak- field limit. Our formulation is based on the standard 
York-Lichnerowicz conformal decomposition and split of the extrinsic curvature into 
its trace and trace-free parts. In this sense it is similar in spirit, if not in detail, to the 
standard approach to the initial value problem [H| and several modern formulations of 
the evolution equations |H1 EI ■ 

The underlying motivation for the new formulation presented in this paper is the 
belief that it is the violation of the constraint equations, and resultant generation 
of spurious energy-momentum sources, which lead to instabilities in numerical 
implementations of Einstein's equations. By solving the Hamiltonian and momentum 
constraints directly, especially in regions of the domain where the gravitational fields 
are strongly dynamic, we automatically prevent the run-away errors which typically 
appear in the Hamiltonian constraint. By automatically satisfying the constraints, we 
guarantee conservation of energy-momentum. 

The formulation of the constraints described in this paper is closest to that 
of Bona, Masso et al [31 01 |S1 , who use the momentum constraints as evolution 
equations. In their early papers the constraints were used purely to eliminate 
spatial derivatives of the extrinsic curvature, thereby casting the equations in strongly 
hyperbolic form Only later are the momentum constraints explicitly described as 
evolution equations [H], although they had always played that role in the formulation. 
More recently Bona et al |2] have incorporated the constraints into the evolution 
system by expanding the Einstein equations themselves, thereby maintaining general 
covariance. Our current approach follows the spirit of their earlier, non-covariant 
formulations. 

We show that the BSSN formulation [HI E3 shares important features with the 
Bona-Masso (BM) approach, in particular their common treatment of the momentum 
constraints. We demonstrate that the BSSN equations incorporate the momentum 
constraints as evolution equations, and we highlight similarities in the way the 
momentum constraints are treated in the BM and BSSN approaches. We then develop 
our own formulation of the momentum constraints, as well as a new approach to 
rewriting the Hamiltonian constraint as an evolution equation. The goal of this 
formulation is the automatic conservation of energy-momentum which we achieve by 
directly solving the constraint equations, thus guaranteeing that they are satisfied 
throughout the evolution. 

We proceed as follows. In the next section we briefly discuss the constrained 
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nature of general relativity, followed by an outline of the approach taken by Bona 
and Masso in treating the momentum constraints. Turning to the BSSN approach, 
we highlight the momentum-conserving properties of the algorithm, before developing 
our own approach. We rewrite both the Hamiltonian and momentum constraints 
as evolution equations, and demonstrate how transparently these equations can 
be incorporated into existing evolution algorithms. Finally, we consider both the 
advantages and the issues raised by our results. 

2. Constraints in General Relativity 

The constraint equations in Einstein's 1915 geometric theory of gravitation play a 
cornerstone role - general relativity is a fully constrained theory. If the four constraints 
are satisfied at every point over every possible spacelike hypersurface of a spacetime, 
then, necessarily, the entire spacetime is a solution of all Einstein's equations. It 
is surprising that the constraints have not played a more pivotal role in numerical 
relativity, given their central importance in theoretical developments of the theory. 
The four constraint equations per point in spacetime are 



known as the Hamiltonian and momentum constraints, respectively. The extrinsic 
curvature is defined as 



with p and Si representing the energy-momentum source terms. In the remainder of 
this paper we work in vacuum (p = S'i = 0), although our analysis is equally valid in 
the presence of matter. 

As it stands, these constraints relate the six components each of gij and Kij, 
and represent an initial value problem which must be solved to obtain consistent 
initial data. Mathematically the constraints are conserved by the evolution equations, 
implying that if they are satisfied on one slice of a foliation, they will be satisfied on all 
future slices. However, it is well know that computational implementations suffer from 
serious errors which propagate in one or more of the constraints. It is this problem 
which we hope to overcome by recasting the constraints as evolution equations. 

3. Bona-Masso treatment of the momentum constraints 

It is our goal to recast all four constraint equations as evolution equations by 
introducing a set of new variables. The goal is to expand the system of evolution 
equations to include the traditional constraints, which necessarily introduces a set 
of new "subsidiary constraints" which must be monitored during evolution. These 
subsidiary constraints take the place of the original Hamiltonian and momentum 
constraints. 

Such an approach has already been taken, perhaps indirectly, by Bona and Masso 
et al while developing hyperbolic formulations of the evolution equations S El El ■ 
Using the momentum constraints to ensure the hyperbolicity of the full set of first- 
order, flux-conservative evolution equations |4|, the Bona-Masso (BM) formulation 



H =R+ {TrKf - K,jK'^ - 2p = 



(1) 
(2) 



(3) 
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introduces a new set of variables 

V^^lg"'{^.gJk~^,g.k). (4) 

Evolution equations for the Vi can be obtained by differentiation and simplification 
using the gij equation together with the momentum constraints 4 . However, the same 
result can be obtained directly from the momentum constraints themselves, which take 
the form jGj: 

dtVk = ~2K'^ {d,,k - S,kd\^) + K'' {dk^J - S,kd^/) (5) 
where we have taken = 1, A^* = for brevity, and 

dkij = ^ dkgtj (6) 

are the first-order flux variables required to reduce the spatial order of the Einstein 
equations. In this approach the definition of Vi, equation Q), is relegated to the 
status of a constraint relating the independent variables Vk, gij and dkij- These are 
the subsidiary constraints in the Bona-Masso approach. 

Bona and Masso provide an elegant incorporation of the momentum constraints 
into the set of dynamical equations, although the four new subsidiary constraints 
are still of a moderately complex form, with no obvious manner of enforcing them 
naturally within an evolution scheme. Before describing an approach to incorporating 
the Hamiltonian constraint into the evolution equations, we pause to consider the 
status of the constraint equations in the BSSN formulation. 



4. Momentum conservation in the BSSN formulation 

As outlined above. Bona and Masso have shown that the momentum constraints can be 
rewritten as first-order evolution equations, which take the form ((SJ when A^ = 1 and 
= 0. They introduce the new variable Vi, defined by equation which appears 
naturally in the momentum constraints when one commutes the spatial derivatives of 
Kij with the time derivatives of gij implicit in the definition of extrinsic curvature. 
This at once removes troublesome spatial derivatives, and explains why eliminating 
spatial derivatives in the evolution equation for Vi with the momentum constraints is in 
fact equivalent to using the momentum constraints themselves as evolution equations. 
In apparent contrast, the BSSN formulation introduces the new variables 

= , (7) 

where g^^ is the contravariant conformal metric. An evolution equation is then 
obtained by differentiation and commutation of spatial and temporal derivatives. 
However, the BSSN formulation uses precisely the same approach to eliminating 
spatial derivatives of the conformal, trace-free portion of the extrinsic curvature A^^ 
as described by Bona and Masso 0. In doing so, the BSSN evolution equation for 
becomes equivalent to the evolution equation obtained directly from the momentum 
constraints themselves, in precisely the same manner as the BM formulation. 

In fact, the BM variable Vi and the BSSN are closely related. Under the 
conformal decomposition 

gij = e^'^'gij with det{gij) = 1 (8) 

the BM variables become 



V^V+4d,(j) 



(9) 
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where Vi is defined like Vi, but in terms of tlie confornial metric g^ . We see tliat Vi 
splits naturally into a portion Vk determined entirely by the conformal geometry, and 
the remainder which depends on the scale factor. A similar relationship was obtained 
previously in the case of a static conformal factor 5 . 

Expanding the conformal portion of Vi , and using the constraint det g = I, we 
find that 

v, = -l~9^Jr\ (10) 

clearly showing that the BSSN variable F* is just the conformal part of the BM variable 
y*. With this realization, it is straightforward to use the momentum constraints, in 
the form of equation Q, to obtain an evolution equation for F*. Not surprisingly, 
the resulting equation is precisely the one used in the standard BSSN formulation to 
evolve F*. 

In their extensive analysis and review of existing formulations of the Einstein 
equations, Shinkai and Yoneda note that the advantages of the BSSN system over the 
standard ADM formulation are due entirely to the introduction of the F' variables, and 
the subsequent elimination of spatial derivatives using the momentum constraint . 
In other words, the sole advantage of the BSSN approach is the use of the momentum 
constraints as evolution equations. The BSSN formulation uses the momentum 
constraints to evolve the conformal portion of Vi, and introduces the constraints 
in their place. 

The relative success of the BSSN formulation provides strong motivation for 
incorporating the Hamiltonian constraint into the set of evolution equations. We 
do this below, as well as proposing an alternative formulation of the momentum 
constraints, with the advantage of a set of extremely natural, and very simple, 
subsidiary constraints. As we shall see, our approach to the momentum constraints is 
opposite to the BSSN choice, since we evolve that portion of Vi arising directly from 
the conformal factor. 



5. Energy conservation: the Hamiltonian constraint as evolution equation 

The key to rewriting the Hamiltonian constraint as an evolution equation is performing 
a conformal decomposition on the three-metric. Rewriting the extrinsic curvature in 
terms of its trace and trace-free parts will allow us to use the Hamiltonian constraint 
to evolve the scale factor. 

We begin with a conformal decomposition of the three-metric, 

9^l^e^*g,, (11) 

with det gij = 1, and split the extrinsic curvature into its trace and trace-free parts, 

A'y = Ay + i.gyTr/^ (12) 

where Tt A = 0. This allows us to rewrite the Hamiltonian constraint as 

H = R.- A,A'^ + I (Tr - 2p, (13) 

where R is the full three-dimensional Ricci scalar calculated from the physical three- 
metric gij . It is therefore a function of (j) and gij , together with their spatial derivatives 
It is not a function of 0, and thus does not play a vital role in the current analysis. 
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We could in principle expand R, as York does, into terms containing spatial derivatives 
of (j) together with the conformal curvature R{g). 

Our aim is to obtain an evolution equation for the conformal factor from the 
Hamiltonian constraint, which requires an understanding of where (j) appears in the 
constraint. Writing 

A., = (^SfS'^~l9,,9^'^Kab (14) 

and applying the conformal decomposition to the definition of extrinsic curvature, 
equation (O, we have 

Kab = [e^'^dtgab - 2V(a7V,) + Agab<p) (15) 

and since 

StS'^~lg^Jg''''^gab^o, (16) 

it is clear that Aij does not depend directly on the time development of (p. The only 
functional dependence on in the Hamiltonian constraint is thus within the Tr K 
term. 

Under the conformal decomposition l|ll|l the trace of the extrinsic curvature 
becomes 

Trif^lv.TV^-^, (17) 

and noting that the Lie derivative of the conformal factor along the shift vector iV* is 
given by 

C^cb = N'dk^ + I dkN^ = iVfciV^ (18) 
6 6 

the trace of the extrinsic curvature can be written as 

TrK^-^{dt-CN)<i^. (19) 

This equation is often used to directly evolve the conformal factor ^U]. Instead, we 
treat it as a definition of Tr K in terms of the time development of the conformal 
factor, allowing us to rewrite the Hamiltonian constraint as the evolution equation for 

The time development of the conformal factor only enters the Hamiltonian 
constraint through the term quadratic in Tr/sT. Proceeding formally, we use equation 
H19|l to rewrite the Hamiltonian constraint ?i = as 



[dt - rjv) = ±y y ^ {A,,A^^ -R + 2p). (20) 

This is a problematic result. The square root causes strife whenever those terms 
within it approach zero (potential computational problems) or become negative. In 
computational tests we find that errors which typically appear as violations of the 
Hamiltonian constraint when using the ADM evolution equations do indeed reappear 
as problems "under the square root" . 

There are several ways of proceeding from the Hamiltonian constraint H13() . In 
general there is no need to replace both Tr K terms, since we wish to maintain linearity 
in the time derivative of the conformal factor. We can therefore write 
(P, r \^ ^ ( R{g)-A,,A^^-2p \ 

{dt-CN)4> = ^^ ^ J, (21) 
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providing an evolution equation for the confornial factor which does not suffer from 
the problems mentioned above. However, problems can still arise if the denominator 
is zero. 

We note that equation H21|l can be rewritten as an evolution equation for a new 
variable 

(-) 

with the advantage that the evolution equation itself is singularity free when Tr K ^ 0. 
Although this form may be more advantageous in computations, it merely relocates 
the problem to the calculation of once ^ is known. 

Another alternative is to monitor Tr K within the code, replacing equation 
H21|l with an alternative whenever the absolute magnitude of Tr K falls below some 
threshold. For example, equation l|19() implies that 

{dt-CN)(b^O (23) 

whenever the trace of the extrinsic curvature is zero. For the remainder of this paper 
we will use equation H21() . 



6. Momentum conservation 



We now turn to the momentum constraints, rewritten in terms of the conformal 
decomposition described in the previous section. Our aim is to find a new set of 
variables which can be evolved using the momentum constraints, and which have a 
simpler subsidiary constraint structure than the BM formulation, equation Q). 
The momentum constraints are 

n, = 9, A-^- + ri,,K) - T%K\ - d, Tr K (24) 

into which we can substitute the definition of AT'^- expressed in terms of its trace and 
trace-free parts. Combining the results in the previous section we find that 

= + ^ 9^J i^kN'^ - 6,^) , (25) 
and thus the momentum constraints are 

The plan is to define the new variable 

= d,c^, (26) 

and commute the temporal and spatial partial derivatives in the momentum constraint 
to obtain an evolution equation for . We note that $j can be viewed as that portion 
of Bona-Masso's Vi variable which is derived purely from the scale factor, as shown 
by equation ©. This is opposite to the choice made in the BSSN formulation, where 
the momentum constraints are used to evolve the conformal portion of Vi. 

Continuing in the manner outlined above, the momentum constraints take the 

form 

dt^j = \do^kN'' -\^^K d,N - ^ V,A^^-, (27) 
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which provides an evolution equation for ^j. However, this can be recast in the more 
convenient form 

{dt - Cn) *j = IdjdkN'' - ^ Tr A' d,N ~ ^ V,^'^ , (28) 

by expanding the covariant derivative of the shift. This is the desired evohition 
equation, derived from the momentum constraints, for the new variables ^j. 

The stability properties of a formulation involving this evolution equation for $j 
remain to be investigated. The BSSN formulation explicitly removes spatial derivatives 
of A^^ (or its conformal part) to improve stability. However, we argue that it is the 
fact that this procedure introduces the momentum constraint as an evolution equation 
which improves the overall stability, not simply the removal of spatial derivatives. 



7. Implementing the "constrained evolution" algorithm 

In the previous sections we described how the standard constraint equations can be 
recast as evolution equations for the conformal factor and its spatial derivatives. It is 
not hard to see that these new forms of the constraints can be easily "bolted onto" 
existing formulations of the evolution equations. 

We first consider the classic g-K, or ADM, formulation of the Einstein equations 
[T^. Incorporating equations and (PSl) into the ADM system, 

i^t-CN)9^J = ~2NK,j (29) 
{dt - Cn) = N - 2KuK'^ + Tr KK.^) - V, V,iV, (30) 

can be approached in a number of ways. We could, for example, introduce Aij, TrK, 
(j) and $j as auxiliary variables, continuing to treat the physical metric and extrinsic 
curvature as the fundamental variables. This is not the most straightforward approach, 
and in general, it would seem advantageous to explicitly construct a conformal, trace- 
free generalization of the ADM equations. 

This approach evolves <j) and gij in place of gij , and replaces Kij with Tr K and 
Aij. The resulting equations are identical to the equivalent BSSN equations, namely 

{dt^CN)~g.j = -2NA,j (31) 

[dt ~ Cn) Ay = e-4<^ [NR., ~ V,V,7V]^^ (32) 

+ N (tt K - 2A,kA\^ (33) 

where Aij = e^*'^Aij and "TF" denotes the trace-free portion of the bracketed 
expression. Various expressions can be derived to evolve Tr K, including the standard 
BSSN equation 

1 .,,2 



[dt- Cn)TtK = N |AyA^ -f - (Tri^)"| - V.V'iV + 2Np. (34) 

The conformal factor (f> and its spatial derivatives ^j are then evolved using the 
Hamiltonian and momentum constraints, equations (|21ll and H28|) : 

N (R{g)-A,jAi -2p 



id.~CN)c^ (35) 

[dt ~ Cn) *j = \d,dkN^ - i TrKd.N - jV,A^. (36) 
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This formulation differs from the existing BSSN approach in two ways. First, 
the momentum "constraint" is used to evolve the variables $j in place of F*. Second, 
the conformal factor is evolved using the Hamiltonian "constraint" , rather than being 
driven directly by the trace of the extrinsic curvature through equation H19|) . 

Similarly, an evolution algorithm based more directly on the BSSN formulation 
can be obtained. The standard BSSN evolution equations are used to evolve gij, 
Aij and TrK, as described above, together with the conformal connection functions 
F' |1(J| . Since these conformal connection functions are already evolved using the 
momentum constraints in the BSSN formulation, our goal of using all four "constraint 
equations" as evolution equations can be achieved by incorporating the Hamiltonian 
equations into the evolution scheme with minor modifications. 

Our form of the Hamiltonian constraint, equation (|21|) . is used to replace the 
standard BSSN evolution equation for the conformal factor which is simply equation 
(|19ll . In this extended-BSSN approach the fundamental variables are still the original 
set, namely (j),gij,Aij,TT K and F*, with the only alteration being that the conformal 
factor is evolved using (|21|l rather than H19|l . This fully-constrained approach is a 
trivial addition to the standard BSSN formulation. 

8. Discussion 

By recasting the standard Hamiltonian and momentum constraints as evolution 
equations, we have guaranteed that the evolution scheme conserves energy- momentum. 
We have also introduced the four subsidiary constraint equations 

det g^j = 1 (37) 

= (38) 

in place of the original Hamiltonian and momentum constraints. These subsidiary 
constraints, or consistency conditions, ensure that the derived variables <i>j and 4>, 
which are evolved independently, continue their expected relation to the original set 
of physical variables. 

The new subsidiary constraints are much simpler than the original Hamiltonian 
and momentum constraints. Codes which violate the Hamiltonian and momentum 
constraints introduce, through computational errors, spurious energy-momentum 
sources. However, violation of the subsidiary constraints does not result in the 
same physical problem, since the original constraints (now in the form of evolution 
equations) are still satisfied to machine precision. Thus no spurious energy-momentum 
can be injected into the system by purely computational errors. It seems likely that 
this advantage explains much of the success of the BSSN formulation, in comparison 
with the ADM form of the equations. In this sense, the BSSN approach may be 
viewed as a natural generalization of the ADM evolution equations to incorporate the 
momentum constraints. 

The major difficulty with the constraint-evolution equations we have obtained is 
the potential for singularities to develop in the Hamiltonian, arising from zeros in the 
denominator. Although this problem can be formally avoided by writing the equation 
in terms of a new variable, the conformal factor must still be calculated from this 
new variable. As such, this approach merely removes singularities from the evolution 
equations without wholly avoiding the problem. One solution, at least for black hole 
spacetimes, would be to use constant-crunch {TrK = constant ^ 0) slicings, which 
have been shown to provide potentially stable foliations of black hole spacetimes ^31 ■ 
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This is, however, not an entirely satisfactory solution, and we continue to investigate 
this issue. 

Work is currently under way to investigate the relative stability of this new 
formulation of the Einstein equations. In particular, we are interested in the relative 
merits of the various forms of the momentum equations (BSSN BM Vi, We 
are also exploring the integration of the full set of constraints into the Bona-Masso 
and BSSN formalisms to determine their mathematical structure, and performing 
computational tests of the new constraint equations. 
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